graph TD
xr["$$\Large \mathbf{x}^{(r)}_{t, i}$$"] --> lambda1(("$$\Large \lambda_{1, t, i}$$"));
betar(("$$\Large \boldsymbol{\beta}_r$$")) --> lambda1;
zr(("$$\Large z^{(r)}_{t, i}$$")) --> lambda1;
drmr: A Bayesian approach to Dynamic Range Models in R
Available at lcgodoy.me/slides/2025-esa/
2025-07-29
Critical Challenge: Predicting species’ responses to global environmental change is vital for conservation and management.
Usual Tool: Species Distribution Models (SDMs) have been the workhorse, correlating occurrences with environmental variables.
Limitations of SDMs:
DRMs explicitly incorporate demographic processes that drive range dynamics (Pagel and Schurr 2012).
Allows linking environmental drivers directly to specific processes.
Potential for more robust forecasting under novel conditions.
Why are they rarely in practice? Despite conceptual appeal, DRMs have been underutilized due to their complexity and computational challenges to fit these models.
drmr: Making DRMs AccessibleGoal: Bridge the gap between DRM potential and practical application.
Key Features:
Stan via cmdstanr for efficient fitting (Gabry et al. 2024).\(Y_{a, t, i}\): Unobserved density of individuals of age \(a\), time \(t\) and site \(i\).
\(Y_{t, i} = \sum_{a} Y_{a, t, i}\): Observed density of individuals (all ages) at time \(t\) and site \(i\).
\(\lambda_{a, t, i}\): Expected age-specific density.
\(\mu_{t, i} = \sum_{a} \lambda_{a, t, i}\): Expected total density.
\(\rho_{t, i} = \mathrm{P}(Y_{t, i} = 0)\): Probability of absence.
\[ f(y_{t, i} \mid \mu_{t, i}, \phi, \rho_{t, i}) = \begin{cases} \rho_{t, i}, & \text{ if } y_{t, i} = 0, \\ (1 - \rho_{t, i}) g(y_{t, i} \mid \mu_{t, i}, \phi), & \text{ if } y_{t, i} > 0. \end{cases} \]
\(y_{t, i}\) is a realization of \(Y_{t, i}\).
\(g(\cdot \mid \mu_{t, i}, \phi)\) is the pdf of a continuous probability distribution.
\(\phi\) is a nuisance parameter.
graph TD
xr["$$\Large \mathbf{x}^{(r)}_{t, i}$$"] --> lambda1(("$$\Large \lambda_{1, t, i}$$"));
betar(("$$\Large \boldsymbol{\beta}_r$$")) --> lambda1;
zr(("$$\Large z^{(r)}_{t, i}$$")) --> lambda1;
graph TD
subgraph "Time t-1"
xs["$$\Large \mathbf{x}^{(s)}_{t - 1, i}$$"] --> s_prev;
betas(("$$\Large \boldsymbol{\beta}_s$$")) --> s_prev;
zs(("$$\Large z^{(s)}_{t - 1, i}$$")) --> s_prev;
f["$$\Large f_{a - 1, t - 1}$$"] --> s_prev;
lambda_prev(("$$\Large \lambda_{a - 1, t - 1, i}$$"))
s_prev(("$$\Large s_{a - 1, t - 1, i}$$"))
end
lambda_prev --> lambda_a(("$$\Large \lambda_{a, t, i}$$"));
s_prev --> lambda_a;
graph TD
%% Mean calculation
lambda_a(("$$\Large \lambda_{a, t, i}$$"));
lambda_a -- $$\Large \sum_a \lambda_{a, t, i}$$ --> mu(("$$\Large \mu_{t, i}$$"));
subgraph "P(Y = 0)"
xt["$$\Large \mathbf{x}^{(p)}_{t, i}$$"] --> rho(("$$\Large \rho_{t, i}$$"));
betat(("$$\Large \boldsymbol{\beta}_p$$")) --> rho;
end
%% Final Observation
mu --> Y{"$$\Large Y_{t, i}$$"};
%% Other components
phi(("$$\Large \phi$$")) --> Y;
rho --> Y;
graph TD
subgraph "Demographic processes"
%% Inputs and Parameters for lambda1
xr["$$\mathbf{x}^{(r)}_{t, i}$$"] --> lambda1(("$$\lambda_{1, t, i}$$"));
betar(("$$\boldsymbol{\beta}_r$$")) --> lambda1;
zr(("$$z^{(r)}_{t, i}$$")) --> lambda1;
%% Recursive part for lambda_a
subgraph "Time t-1"
xs["$$\mathbf{x}^{(s)}_{t - 1, i}$$"] --> s_prev;
betas(("$$\boldsymbol{\beta}_s$$")) --> s_prev;
zs(("$$z^{(s)}_{t - 1, i}$$")) --> s_prev;
f["$$f_{a - 1, t - 1}$$"] --> s_prev;
lambda_prev(("$$\lambda_{a - 1, t - 1, i}$$"))
s_prev(("$$s_{a - 1, t - 1, i}$$"))
end
lambda_prev --> lambda_a(("$$\lambda_{a, t, i}$$"));
s_prev --> lambda_a;
end
%% Mean calculation
lambda_a --> mu(("$$\mu_{t, i}$$"));
subgraph "Zero-augmentation"
xt["$$\mathbf{x}^{(p)}_{t, i}$$"] --> rho(("$$\rho_{t, i}$$"));
betat(("$$\boldsymbol{\beta}_p$$")) --> rho;
end
%% Final Observation
mu --> Y{"$$Y_{t, i}$$"};
%% Other components
phi(("$$\phi$$")) --> Y;
rho --> Y;
An example analysis uses Summer flounder (Paralichthys dentatus) data from 1982-2016 NOAA bottom trawl surveys (Fredston et al. 2025) to illustrate the package’s features.
The data spans the US Atlantic coast (Cape Hatteras, NC to the Canada/Maine border) and was aggregated from individual hauls into 10 latitudinal patches with varying areas.
Response variable: Density (count per unit area).
Environmental drivers: SST and SBT.
| Model | Non-neg. PDF | \(\rho_{i, t}\) | Recruitment (or Density) | Survival |
|---|---|---|---|---|
| DRM (rec) | Gamma | Effort + estimated density | SST + SST2 + AR(1) | Intercept + \(f_{a, t}\) |
| DRM (surv) | Gamma | Effort + estimated density | AR(1) | SBT + SBT2 + \(f_{a, t}\) |
| SDM (GLM) | LogNormal | Effort | SST + SST2 |
drm_rec <-
fit_drm(.data = dat_train,
y_col = "dens", ## response variable: density
time_col = "year", ## vector of time points
site_col = "patch",
family = "gamma", ## other options: lognormal, loglogistic
seed = 202505,
formula_zero = ~ 1 + c_hauls,
formula_rec = ~ 1 + c_sst + I(c_sst * c_sst),
formula_surv = ~ 1,
f_mort = f_train,
n_ages = NROW(f_train),
adj_mat = adj_mat, ## A matrix for movement routine
ages_movement = c(0, 0,
rep(1, 12),
0, 0), ## ages allowed to move
.toggles = list(ar_re = "rec", ## other options: "surv", "dens"
movement = 1,
est_surv = 1,
est_init = 0,
minit = 1),
.priors = list(pr_phi_a = 1, pr_phi_b = .1,
pr_alpha_a = 4.2, pr_alpha_b = 5.8,
pr_zeta_a = 7, pr_zeta_b = 3))drm_surv <-
fit_drm(.data = dat_train,
y_col = "dens", ## response variable: density
time_col = "year", ## vector of time points
site_col = "patch",
family = "gamma", ## other options: lognormal, loglogistic
seed = 202505,
formula_zero = ~ 1 + c_hauls,
formula_rec = ~ 1,
formula_surv = ~ 1 + c_sbt + I(c_sst * c_sbt),
f_mort = f_train,
n_ages = NROW(f_train),
adj_mat = adj_mat, ## A matrix for movement routine
ages_movement = c(0, 0,
rep(1, 12),
0, 0), ## ages allowed to move
.toggles = list(ar_re = "rec", ## other options: "surv", "dens"
movement = 1,
est_surv = 1,
est_init = 0,
minit = 1),
.priors = list(pr_phi_a = 1, pr_phi_b = .1,
pr_alpha_a = 4.2, pr_alpha_b = 5.8,
pr_zeta_a = 7, pr_zeta_b = 3))predict_drm function requires at least three arguments:
drm: Output of a fit_drm callnew_data: A data.frame where we seek to obtain predictionsseed: For reproducibilityWhen \(f_{a, t}\) is used at the time of model fit, we also need to provide f_test to predict_drm.
past_data is necessary when survival is a function of environmental variables.
The drmr substantially lowers the barrier for ecologists to use the DRM in their applications.
The code is easy to use and takes advantage of what has been developed for Stan: visualization, diagnostic tools, and estimation.
drmr allows for empirically testing which processes are more important to predict the distribution of a species.
The more complex a model is, the more (and better) data we need to be able to estimate those relationships.
Include population dynamic models that explicitly relate adult density to recruitment (e.g., Ricker, Belverton-Holt)
GAM-like non-parametric relationships between processes and the environment
Support for length-composition data
More realistic movement routines